// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                           avtHistogramFilter.C                            //
// ************************************************************************* //

#include <avtHistogramFilter.h>

#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkPointDataToCellData.h>
#include <vtkPolyData.h>
#include <vtkUnsignedIntArray.h>

#include <PickAttributes.h>
#include <PickVarInfo.h>

#include <avtDataAttributes.h>
#include <avtDatasetExaminer.h>
#include <avtExtents.h>
#include <avtParallel.h>
#include <avtOriginatingSource.h>

#include <BadCellException.h>
#include <ImproperUseException.h>

#include <float.h>

#include <string>
#include <vector>
#include <avtExecutionManager.h>

using std::string;
using std::vector;

// ****************************************************************************
//  Method: avtHistogramFilter constructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Jun 26 09:04:54 PDT 2003
//
// ****************************************************************************

avtHistogramFilter::avtHistogramFilter()
{
    bins = NULL;
}


// ****************************************************************************
//  Method: avtHistogramFilter destructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Jun 26 09:04:54 PDT 2003
//
// ****************************************************************************

avtHistogramFilter::~avtHistogramFilter()
{
    if (bins != NULL)
        delete [] bins;

    VisitMutexDestroy("avtHistogramFilter");
}


// ****************************************************************************
//  Method: avtHistogramFilter::SetAttributes
//
//  Purpose:
//      Sets the attributes for this filter.
//
//  Programmer: Hank Childs
//  Creation:   June 26, 2003
//
// ****************************************************************************

void
avtHistogramFilter::SetAttributes(const HistogramAttributes &h_atts)
{
    atts = h_atts;
}


// ****************************************************************************
//  Method: avtHistogramFilter::PreExecute
//
//  Purpose:
//      Initializes for our actual execute.  Sets up the bins, etc.
//
//  Programmer: Hank Childs
//  Creation:   June 26, 2003
//
//  Modifications:
//
//    Hank Childs, Wed May 24 10:18:09 PDT 2006
//    Add support for array variables.
//    
//    Dave Pugmire, Thu Nov 01 12:39:07 EDT 2007
//    Support for log, sqrt scaling.
//
//    Hank Childs, Wed Mar  5 09:54:20 PST 2008
//    Don't get the data range if the min and max are already specified.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::PreExecute(void)
{
    avtDataTreeIterator::PreExecute();

    InputSetActiveVariable(pipelineVariable);

    if (atts.GetBasedOn() == HistogramAttributes::ManyZonesForSingleVar)
    {
        if ( !atts.GetMinFlag() || !atts.GetMaxFlag() )
        {
            // Original data extents can be accessed directly via GetDataExtents.
            // ActualDataExtents may not be set by some operators, so we use the
            // DataExaminer to get the actual extents when the user selects
            // "Current Plot".
            if (atts.GetLimitsMode() == HistogramAttributes::LimitsMode::OriginalData)
            {
                GetDataExtents(dataValueRange, pipelineVariable);
            }
            else if (atts.GetLimitsMode() == HistogramAttributes::LimitsMode::CurrentPlot)
            {
                avtDataset_p ds = GetTypedInput();
                avtDatasetExaminer::GetDataExtents(ds, dataValueRange);
            }
        }

        SetWorkingMin( (atts.GetMinFlag() ? atts.GetMin() : dataValueRange[0]) );
        SetWorkingMax( (atts.GetMaxFlag() ? atts.GetMax() : dataValueRange[1]) );
        SetWorkingNumBins( atts.GetNumBins() );
        
        if (bins != NULL)
            delete [] bins;
        bins = new double[workingNumBins];
        for (int i = 0 ; i < workingNumBins ; i++)
            bins[i] = 0.;
    }
    else
        workingNumBins = 0;
}


// ****************************************************************************
//  Method: avtHistogramFilter::PostExecute
//
//  Purpose:
//      Sets up the final curve based on the bins.
//
//  Programmer: Hank Childs
//  Creation:   June 26, 2003
//
//  Modifications:
//
//    Hank Childs, Sat Oct 18 11:14:06 PDT 2003
//    Make sure the limits of the histogram plot always have y=0.
//
//    Hank Childs, Sat Oct 18 11:45:13 PDT 2003
//    Make sure that the bins are of a uniform size.
//
//    Brad Whitlock, Thu Jul 22 17:02:11 PST 2004
//    Changed the xlabel and xunits in the output data attributes so the
//    histogram plot would have the units along the axes.
//
//    Kathleen Bonnell, Thu Jan  6 10:34:57 PST 2005 
//    Remove TRY-CATCH block in favor of testing for ValidActiveVariable.
//
//    Kathleen Bonnell, Fri May 13 11:10:35 PDT 2005 
//    Fix memory leak. 
//
//    Hank Childs, Wed May 24 11:27:24 PDT 2006
//    Add support for array variables.
//
//    Hank Childs, Thu Sep 14 09:16:23 PDT 2006
//    Fix indexing bug pointed out by Matt Wheeler.
//
//    Hank Childs, Fri Jan  5 11:49:31 PST 2007
//    Reverse indexing for bins that have negative values.
//
//    Hank Childs, Fri Jan 12 15:27:01 PST 2007
//    Use bin spacing if possible.
//
//    Cyrus Harrison, Wed Mar  7 15:46:29 PST 2007
//    Added support for point histograms and true "Frequency" histograms
//
//    Jeremy Meredith, Wed Mar 14 11:14:56 EDT 2007
//    Call avtDataTreeIterator::PostExecute instead of PreExecute.
//
//    Dave Pugmire, Thu Nov 01 12:39:07 EDT 2007
//    Support for log, sqrt scaling.    
//
//    Hank Childs, Wed Dec 12 21:59:49 PST 2007
//    Add support for weighting by a variable.
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
//    Kevin Bensema, Mon Nov 11 13:04 PST 2013
//    Added code to change the y-axis units to "Probability" if the histogram
//    has been normalized.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::PostExecute(void)
{
    int  i;

    avtDataTreeIterator::PostExecute();

    if (atts.GetBasedOn() == HistogramAttributes::ManyVarsForSingleZone)
    {
        //
        // Initialize for all the processors that got no data.
        //
        workingNumBins = UnifyMaximumValue(workingNumBins);
        if (workingNumBins == 0)
        {
            EXCEPTION1(VisItException, "The histogram could not be generated."
                             "  The most likely source of error is that the "
                             "zone id or domain number is invalid.  If this"
                             " is not the case, please contact a VisIt"
                             " developer.");
        }
        if (bins == NULL)
        {
            bins = new double[workingNumBins];
            for (i = 0 ; i < workingNumBins ; i++)
                bins[i] = 0.;
        }
        workingRange[0] = 0.;
        workingRange[1] = (double)workingNumBins;
    }

    double *newBins = new double[workingNumBins];
    SumDoubleArrayAcrossAllProcessors(bins, newBins, workingNumBins);
    for (i = 0 ; i < workingNumBins ; i++)
        bins[i] = newBins[i];

    ScaleBins();

    delete [] newBins;

    if (PAR_Rank() != 0)
    {
        SetOutputDataTree(new avtDataTree());
        return;
    }

    vtkPolyData *output = vtkPolyData::New();

    bool spaceBins = false;
    vector<double> ranges;
    if (atts.GetBasedOn() == atts.ManyVarsForSingleZone)
        if (atts.GetUseBinWidths() == true)
        {
            ranges = 
                 GetInput()->GetInfo().GetAttributes().GetVariableBinRanges();
            if (ranges.size() == (size_t)workingNumBins+1)
                spaceBins = true;
        }

    if (atts.GetOutputType() == HistogramAttributes::Curve)
    {
        vtkPoints *pts = vtkPoints::New();
        pts->SetNumberOfPoints(workingNumBins);
        double jump = (GetWorkingMax() - GetWorkingMin()) / (workingNumBins);
        for (i = 0 ; i < workingNumBins ; i++)
        {
            double pt[3];
            if (spaceBins)
                pt[0] = (ranges[i]+ranges[i+1])/2.;
            else
                pt[0] = GetWorkingMin() + (i+0.5)*jump;
            pt[1] = bins[i];
            pt[2] = 0.;
            pts->SetPoint(i, pt);
        }
    
        output->SetPoints(pts);
        pts->Delete();
    
        int nSegments = workingNumBins-1;
        output->Allocate(nSegments);
        vtkIdType line[2];
        for (int i = 0 ; i < nSegments ; i++)
        {
            line[0] = i;
            line[1] = i+1;
            output->InsertNextCell(VTK_LINE, 2, line);
        }
    }
    else
    {
        vtkPoints *pts = vtkPoints::New();
        pts->SetNumberOfPoints(4*workingNumBins);
        double jump = (GetWorkingMax() - GetWorkingMin()) / (workingNumBins);

        int ptIndex = 0;
        double pt[3] = { 0., 0., 0. };

        for (i = 0 ; i < workingNumBins ; i++)
        {
            //
            // Do the last two points for bin i.
            //
            if (spaceBins)
                pt[0] = ranges[i];
            else
                pt[0] = jump*i + GetWorkingMin();
            pt[1] = 0.;
            pts->SetPoint(ptIndex++, pt);

            pt[1] = bins[i];
            pts->SetPoint(ptIndex++, pt);
 
            //
            // Now do the first two points for bin i+1.
            //
            if (spaceBins)
                pt[0] = ranges[i+1];
            else
                pt[0] = jump*(i+1) + GetWorkingMin();
            pt[1] = 0.;
            pts->SetPoint(ptIndex++, pt);

            pt[1] = bins[i];
            pts->SetPoint(ptIndex++, pt);
        }
    
        output->SetPoints(pts);
        pts->Delete();
    
        int ncells = 5*workingNumBins;
        output->Allocate(ncells);

        //
        // Add the line segments for each bin.
        //
        vtkIdType line[2];
        for (i = 0 ; i < workingNumBins ; i++)
        {
            line[0] = 4*i;
            line[1] = 4*i+1;
            output->InsertNextCell(VTK_LINE, 2, line);
            line[0] = 4*i+2;
            line[1] = 4*i+3;
            output->InsertNextCell(VTK_LINE, 2, line);
            line[0] = 4*i;
            line[1] = 4*i+2;
            output->InsertNextCell(VTK_LINE, 2, line);
            line[0] = 4*i+1;
            line[1] = 4*i+3;
            output->InsertNextCell(VTK_LINE, 2, line);
        }

        //
        // Now add the polygons in the middle.
        //
        vtkIdType quad[4];
        for (int i = 0 ; i < workingNumBins ; i++)
        {
            if (bins[i] >= 0.)
            {
                quad[0] = 4*i;
                quad[1] = 4*i+2;
                quad[2] = 4*i+3;
                quad[3] = 4*i+1;
            }
            else
            {
                quad[0] = 4*i;
                quad[1] = 4*i+1;
                quad[2] = 4*i+3;
                quad[3] = 4*i+2;
            }
            output->InsertNextCell(VTK_QUAD, 4, quad);
        }
    }

    SetOutputDataTree(new avtDataTree(output, -1));
    output->Delete();

    //
    // Set up the extents of the output.
    //
    double low = +DBL_MAX;
    double hi  = -DBL_MAX;
    for (i = 0 ; i < workingNumBins ; i++)
    {
        low = (low < bins[i] ? low : bins[i]);
        hi  = (hi  > bins[i] ? hi  : bins[i]);
    }
    double extents[6];
    if (spaceBins)
    {
        extents[0] = ranges[0];
        extents[1] = ranges[workingNumBins];
    }
    else
    {
        extents[0] = GetWorkingMin();
        extents[1] = GetWorkingMax();
    }
    if (atts.GetBasedOn() == HistogramAttributes::ManyZonesForSingleVar)
        extents[2] = 0.;  // We always start from 0.  lo is misleading.
    else
        extents[2] = (low < 0. ? low : 0.);
    extents[3] = hi;
    extents[4] = extents[5] = 0.;

    // Try to do something reasonable if we have constant data.
    if (extents[0] == extents[1])
    {
        extents[0] = GetWorkingMin() - hi/2.;
        extents[1] = GetWorkingMax() + hi/2.;
    }

    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    outAtts.GetOriginalSpatialExtents()->Set(extents);
    outAtts.GetThisProcsOriginalSpatialExtents()->Set(extents);
    outAtts.GetDesiredSpatialExtents()->Set(extents);
    outAtts.GetActualSpatialExtents()->Set(extents);
    outAtts.GetThisProcsActualSpatialExtents()->Set(extents);

    //
    // Set the X-axis's units to match the variable units.
    // 
    if (GetInput()->GetInfo().GetAttributes().ValidActiveVariable())
    {
        if(GetInput()->GetInfo().GetAttributes().GetVariableUnits() != "")
        {
            string xlabel = "";
            if ( atts.GetDataScale() == HistogramAttributes::Linear )
                xlabel = string( "Variable " ) + pipelineVariable;
            else if ( atts.GetDataScale() == HistogramAttributes::Log )
                xlabel = string( "Variable log10(" ) + pipelineVariable + string(") ");
            if ( atts.GetDataScale() == HistogramAttributes::SquareRoot )
                xlabel = string( "Variable sqrt(" ) + pipelineVariable + string(") ");
            
            outAtts.SetXLabel( xlabel );
            outAtts.SetXUnits(GetInput()->GetInfo().GetAttributes().GetVariableUnits());
        }
        else
        {
            string str = "";
            if ( atts.GetDataScale() == HistogramAttributes::Linear )
                str = pipelineVariable;
            else if ( atts.GetDataScale() == HistogramAttributes::Log )
                str = string( "log10(" ) + pipelineVariable + string(") ");
            if ( atts.GetDataScale() == HistogramAttributes::SquareRoot )
                str = string( "sqrt(" ) + pipelineVariable + string(") ");            
            outAtts.SetXUnits(str );
        }
    }

    if (atts.GetBasedOn() == HistogramAttributes::ManyVarsForSingleZone)
    {
        outAtts.SetYLabel("Value");
        outAtts.SetYUnits("");
    }
    else
    {
        int topo = GetInput()->GetInfo().GetAttributes().GetTopologicalDimension();
        string yunits = "";
        if (atts.GetHistogramType() == HistogramAttributes::Frequency)
        {
            if (topo == 0)
                yunits = "# of Points";
            else
                yunits = "# of Cells";

            if(atts.GetNormalizeHistogram() == true)
            {
              yunits = "Probability";
            }
        }
        else if (atts.GetHistogramType() == HistogramAttributes::Weighted)
        {
            // For particles, pretend we have frequency.
            if (topo == 0)
            {
                yunits = "# of Points";
            }
            else if (topo==3)
            {
                yunits = "Volume";
            }
            else // topo = 2 
            {
                if (GetInput()->GetInfo().GetAttributes().GetMeshCoordType() == AVT_XY)
                    yunits = "Area";
                else
                    yunits = "Revolved Volume";
            }
        }
        else
        {
            if (atts.GetWeightVariable() == "default")
                yunits = pipelineVariable;
            else
                yunits = atts.GetWeightVariable();
        }
        string str = yunits;
        if ( atts.GetBinScale() == HistogramAttributes::Log )
            str = "log10(" + yunits + ") ";
        else if ( atts.GetBinScale() == HistogramAttributes::SquareRoot )
            str = "sqrt(" + yunits + ") ";

        outAtts.SetYUnits(str);
    }
}


// ****************************************************************************
//  Method: avtHistogramFilter::ThreadSafe
//
//  Purpose:
//      Returns if this filter is thread safe. There is one mode that we
//  are not thread safe.
//
//  Programmer: David Camp
//  Creation:   March 11 2013
//
// ****************************************************************************
bool
avtHistogramFilter::ThreadSafe(void)
{
    if (atts.GetBasedOn() == HistogramAttributes::ManyVarsForSingleZone)
    {
        return( false );
    }
    return( true );
}

// ****************************************************************************
//  Method: avtHistogramFilter::ExecuteData
//
//  Purpose:
//      Does the actual VTK code to modify the dataset.
//
//  Arguments:
//      inDR      The input data representation.
//
//  Returns:      The output data representation.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Jun 26 09:04:54 PDT 2003
//
//  Modifications:
//
//    Hank Childs, Sat Oct 18 11:41:13 PDT 2003
//    Make the bins be uniform in size.
//
//    Hank Childs, Wed May 24 09:57:40 PDT 2006
//    Add support for taking a histogram based on an array variable.
//  
//    Cyrus Harrison, Thu Mar  8 08:08:38 PST 2007
//    Add support for point histograms and frequency histograms
//
//    Hank Childs, Tue Dec 11 20:25:22 PST 2007
//    Add support for variable weighting.
//
//    Eric Brugger, Tue Aug 19 10:28:44 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtHistogramFilter::ExecuteData(avtDataRepresentation *inDR)
{
    //
    // Get the VTK data set and domain number.
    //
    vtkDataSet *inDS = inDR->GetDataVTK();
    int domain = inDR->GetDomain();

    if (atts.GetBasedOn() == HistogramAttributes::ManyVarsForSingleZone)
    {
        ArrayVarExecute(inDS, domain);
    }
    else if (atts.GetHistogramType() == HistogramAttributes::Frequency)
    {
        FreqzExecute(inDS);
    }
    else if (atts.GetHistogramType() == HistogramAttributes::Weighted)
    {
        if (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() == 0)
            FreqzExecute(inDS);
        else
            WeightedExecute(inDS);
    }
    else if (atts.GetHistogramType() == HistogramAttributes::Variable)
    {
        VariableExecute(inDS);
    }
    return NULL;
}

// ****************************************************************************
//  Method: avtHistogramFilter::FreqzExecute
//
//  Purpose:
//      Standard Frequency Histogram (bins accumlate cell or point counts)
//
//  Notes: Based on Hank's old StandardExecute method
//
//  Programmer: Cyrus Harrison
//  Creation:   March 7, 2007
//
//  Modifications:
//
//    Hank Childs, Mon Oct 22 15:58:59 PDT 2007
//    Ignore ghost data.
//
//    Dave Pugmire, Thu Nov 01 12:39:07 EDT 2007
//    Support for log, sqrt scaling.    
//
//    Hank Childs, Wed Mar  5 10:17:00 PST 2008
//    Don't recenter nodal data to be zonal.
//
//    David Camp, Tue Mar 12 13:19:37 PDT 2013
//    Made function thread safe.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
//
// ****************************************************************************

void
avtHistogramFilter::FreqzExecute(vtkDataSet *inDS)
{
    //
    // Get the variable that we are binning by.
    //
    const char *var = pipelineVariable;
    vtkDataArray *bin_arr = inDS->GetPointData()->GetArray(var);
    if (bin_arr == NULL)
        bin_arr = inDS->GetCellData()->GetArray(var);
    if (bin_arr == NULL)
        EXCEPTION0(ImproperUseException);

    unsigned char *ghosts = NULL;
    if (inDS->GetCellData()->GetArray("avtGhostZones") != NULL)
    {
        vtkUnsignedCharArray *g = (vtkUnsignedCharArray *)
                            inDS->GetCellData()->GetArray("avtGhostZones");
        if (g->GetNumberOfTuples() == bin_arr->GetNumberOfTuples())
            ghosts = g->GetPointer(0);
    }

    double *threadBins;
#if defined(VISIT_THREADS)
    threadBins = new double[workingNumBins];
    memset( threadBins, 0, sizeof(double)*workingNumBins );
#else
    threadBins = bins;
#endif // VISIT_THREADS

    //
    // Now we will walk through each value and sort them into bins.
    //
    int nvals = bin_arr->GetNumberOfTuples();
    for (int i = 0 ; i < nvals ; i++)
    {
        if (ghosts != NULL && ghosts[i] != '\0')
            continue;
        double val = bin_arr->GetTuple1(i);
        int index = ComputeBinIndex( val );
        if ( index < 0 )
            continue;

        if (index >= workingNumBins)
            index = workingNumBins-1;
        threadBins[index]++;
    }

#if defined(VISIT_THREADS)
    VisitMutexLock("avtHistogramFilter");
    for (int i = 0 ; i < workingNumBins ; i++)
    {
        bins[i] += threadBins[i];
    }
    VisitMutexUnlock("avtHistogramFilter");
    delete [] threadBins;
#endif // VISIT_THREADS
}


// ****************************************************************************
//  Method: avtHistogramFilter::WeightedExecute
//
//  Purpose:
//      Weighted Histogram (bins accumlate area of volume of cells)
//
//  Notes: Based on Hank's old StandardExecute method
//
//  Programmer: Cyrus Harrison
//  Creation:   March 7, 2007
//
//  Modifications:
//
//    Hank Childs, Mon Oct 22 15:58:59 PDT 2007
//    Ignore ghost data.
//
//    Dave Pugmire, Thu Nov 01 12:39:07 EDT 2007
//    Support for log, sqrt scaling.    
//
//    David Camp, Tue Mar 12 13:19:37 PDT 2013
//    Made function thread safe.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::WeightedExecute(vtkDataSet *inDS)
{
    // Get the "_amounts".  This is the area or volume that each cell
    // takes up.
    //
    vtkDataArray *amount_arr = inDS->GetCellData()->GetArray("_amounts");
    if (amount_arr == NULL)
        EXCEPTION0(ImproperUseException);
    //
    // Get the variable that we are binning by.
    //
    const char *var = pipelineVariable;
    vtkDataArray *bin_arr = NULL;
    bool ownBinArr = false;

    if(GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() == 0)
    {
        // weighted case does not make sense for point data
        EXCEPTION0(ImproperUseException);
    }
    
    if (inDS->GetPointData()->GetArray(var) != NULL)
    {
        // in the 2d or 3d case make sure to get zone centered data
        //
        // The input is point-centered, but we would prefer zone-centered.
        //
        vtkDataSet *new_in_ds = (vtkDataSet *) inDS->NewInstance();
        new_in_ds->CopyStructure(inDS);
        new_in_ds->GetPointData()->AddArray(
                                          inDS->GetPointData()->GetArray(var));
        vtkPointDataToCellData *pd2cd = vtkPointDataToCellData::New();
        pd2cd->SetInputData(new_in_ds);
        pd2cd->Update();

        bin_arr = pd2cd->GetOutput()->GetCellData()->GetArray(var);
        bin_arr->Register(NULL);
        ownBinArr = true;

        new_in_ds->Delete();
        pd2cd->Delete();
    }
    else
    {
        bin_arr = inDS->GetCellData()->GetArray(var);
    }
    if (bin_arr == NULL)
        EXCEPTION0(ImproperUseException);

    unsigned char *ghosts = NULL;
    if (inDS->GetCellData()->GetArray("avtGhostZones") != NULL)
    {
        vtkUnsignedCharArray *g = (vtkUnsignedCharArray *)
                            inDS->GetCellData()->GetArray("avtGhostZones");
        if (g->GetNumberOfTuples() == bin_arr->GetNumberOfTuples())
            ghosts = g->GetPointer(0);
    }

    //
    // Now we will walk through each value and sort them into bins.
    //
    int nvals = bin_arr->GetNumberOfTuples();
    
    double *threadBins;
#if defined(VISIT_THREADS)
    threadBins = new double[workingNumBins];
    memset( threadBins, 0, sizeof(double)*workingNumBins );
#else
    threadBins = bins;
#endif // VISIT_THREADS

    for (int i = 0 ; i < nvals ; i++)
    {
        if (ghosts != NULL && ghosts[i] != '\0')
            continue;
        double val = bin_arr->GetTuple1(i);
        int index = ComputeBinIndex( val );
        if ( index < 0 )
            continue;
        if (index >= workingNumBins)
            index = workingNumBins-1;
        double amount = amount_arr->GetTuple1(i);
        threadBins[index] += amount;
    }

    if (ownBinArr)
        bin_arr->Delete();

#if defined(VISIT_THREADS)
    VisitMutexLock("avtHistogramFilter");
    for (int i = 0 ; i < workingNumBins ; i++)
    {
        bins[i] += threadBins[i];
    }
    VisitMutexUnlock("avtHistogramFilter");
    delete [] threadBins;
#endif // VISIT_THREADS
}


// ****************************************************************************
//  Method: avtHistogramFilter::VariableExecute
//
//  Purpose:
//      Histogram weighted by a variable
//
//  Programmer: Hank Childs
//  Creation:   December 11, 2007
//
//  Modifications:
//
//    David Camp, Tue Mar 12 13:19:37 PDT 2013
//    Made function thread safe.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::VariableExecute(vtkDataSet *inDS)
{
    const char *histIndexVarName = pipelineVariable;
    const char *weightVarName    = atts.GetWeightVariable().c_str();
    if (strcmp(weightVarName, "default") == 0)
        weightVarName = histIndexVarName;

    bool histIndexVarIsZonal = false;
    bool weightVarIsZonal    = false;

    if (inDS->GetCellData()->GetArray(histIndexVarName) != NULL)
        histIndexVarIsZonal = true;
    else if (inDS->GetPointData()->GetArray(histIndexVarName) == NULL)
    {
        EXCEPTION0(ImproperUseException);
    }
    if (inDS->GetCellData()->GetArray(weightVarName) != NULL)
        weightVarIsZonal = true;
    else if (inDS->GetPointData()->GetArray(weightVarName) == NULL)
    {
        EXCEPTION0(ImproperUseException);
    }

    int           nvals        = 0;
    vtkDataArray *histIndexVar = NULL;
    vtkDataArray *weightVar    = NULL;
    bool          ownWeightVar = false;
    if (histIndexVarIsZonal == weightVarIsZonal)
    {
        if (histIndexVarIsZonal)
        {
            nvals        = inDS->GetNumberOfCells();
            histIndexVar = inDS->GetCellData()->GetArray(histIndexVarName);
            weightVar    = inDS->GetCellData()->GetArray(weightVarName);
        }
        else
        {
            nvals        = inDS->GetNumberOfPoints();
            histIndexVar = inDS->GetPointData()->GetArray(histIndexVarName);
            weightVar    = inDS->GetPointData()->GetArray(weightVarName);
        }
    }
    else
    {
        // Recenter to the the same centering as the histogram index variable.
        vtkDataSet *new_in_ds = (vtkDataSet *) inDS->NewInstance();
        new_in_ds->CopyStructure(inDS);

        if (histIndexVarIsZonal)
        {
            nvals        = inDS->GetNumberOfCells();
            histIndexVar = inDS->GetCellData()->GetArray(histIndexVarName);
            new_in_ds->GetPointData()->AddArray(
                                inDS->GetPointData()->GetArray(weightVarName));
            vtkPointDataToCellData *pd2cd = vtkPointDataToCellData::New();
            pd2cd->SetInputData(new_in_ds);
            pd2cd->Update();
            weightVar = pd2cd->GetOutput()->GetCellData()->GetArray(weightVarName);
            weightVar->Register(NULL);
            pd2cd->Delete();
        }
        else
        {
            nvals        = inDS->GetNumberOfPoints();
            histIndexVar = inDS->GetPointData()->GetArray(histIndexVarName);
            new_in_ds->GetCellData()->AddArray(
                                inDS->GetCellData()->GetArray(weightVarName));
            vtkCellDataToPointData *cd2pd = vtkCellDataToPointData::New();
            cd2pd->SetInputData(new_in_ds);
            cd2pd->Update();
            weightVar = cd2pd->GetOutput()->GetPointData()->GetArray(weightVarName);
            weightVar->Register(NULL);
            cd2pd->Delete();
        }

        new_in_ds->Delete();
        ownWeightVar = true;
    }

    unsigned char *ghosts = NULL;
    if (inDS->GetCellData()->GetArray("avtGhostZones") != NULL)
    {
        vtkUnsignedCharArray *g = (vtkUnsignedCharArray *)
                            inDS->GetCellData()->GetArray("avtGhostZones");
        if (g->GetNumberOfTuples() == nvals)
            ghosts = g->GetPointer(0);
    }

    double *threadBins;
#if defined(VISIT_THREADS)
    threadBins = new double[workingNumBins];
    memset( threadBins, 0, sizeof(double)*workingNumBins );
#else
    threadBins = bins;
#endif // VISIT_THREADS

    //
    // Now we will walk through each value and sort them into bins.
    //
    for (int i = 0 ; i < nvals ; i++)
    {
        if (ghosts != NULL && ghosts[i] != '\0')
            continue;
        double val = histIndexVar->GetTuple1(i);
        int index = ComputeBinIndex( val );
        if ( index < 0 )
            continue;
        if (index >= workingNumBins)
            index = workingNumBins-1;
        double amount = weightVar->GetTuple1(i);
        threadBins[index] += amount;
    }

    if (ownWeightVar)
        weightVar->Delete();

#if defined(VISIT_THREADS)
    VisitMutexLock("avtHistogramFilter");
    for (int i = 0 ; i < workingNumBins ; i++)
    {
        bins[i] += threadBins[i];
    }
    VisitMutexUnlock("avtHistogramFilter");
    delete [] threadBins;
#endif // VISIT_THREADS
}


// ****************************************************************************
//  Method: avtHistogramFilter::ArrayVarExecute
//
//  Purpose:
//      Take a histogram by finding the right zone and creating a histogram
//      from its array variables.
//
//  Programmer: Hank Childs
//  Creation:   May 24, 2006
//
//  Modifications:
//
//    David Camp, Tue Mar 12 13:19:37 PDT 2013
//    Hank and I look at this function. It should be thread safe because only
//    one data block should have the zone in it.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::ArrayVarExecute(vtkDataSet *inDS, int chunk)
{
    int domain = atts.GetDomain();
    int blockOrigin = GetInput()->GetInfo().GetAttributes().GetBlockOrigin();
    domain -= blockOrigin;
    if (chunk != domain)
        return;

    vtkDataArray *arr = inDS->GetCellData()->GetArray(pipelineVariable);
    if (arr == NULL)
    {
        // "Cell" data is stored as point data for point meshes.
        arr = inDS->GetPointData()->GetArray(pipelineVariable);
        if (arr == NULL)
            EXCEPTION0(ImproperUseException);
    }

    vtkUnsignedIntArray *cid = (vtkUnsignedIntArray *) 
                       inDS->GetCellData()->GetArray("avtOriginalCellNumbers");
    if (cid == NULL)
    {
        // This should come down because we requested it in the data
        // specification.
        EXCEPTION0(ImproperUseException);
    }
    int ncomps = cid->GetNumberOfComponents();
    int compOffset = ncomps-1; // 0 if 1 comp, 1 if 2 comps
    unsigned int *ptr = cid->GetPointer(0);

    // Pick returns with cell origin, but avtOriginalCellNumbers does not
    // have a cell origin.  So subtract it off here.
    int zone = atts.GetZone();
    int cellOrigin  = GetInput()->GetInfo().GetAttributes().GetCellOrigin();
    zone -= cellOrigin;
    int ncells = inDS->GetNumberOfCells();
    for (int i = 0 ; i < ncells ; i++)
    {
        if (ptr[ncomps*i + compOffset] == (size_t)zone)
        {
            workingNumBins = arr->GetNumberOfComponents();
            double *vals = new double[workingNumBins];

            arr->GetTuple(i, vals);

            if (bins != NULL)
                delete [] bins;

            bins = new double[workingNumBins];
            for (int i = 0 ; i < workingNumBins ; i++)
            {
                bins[i] = vals[i];
            }
            delete [] vals;
        }
    }
}


// ****************************************************************************
//  Method: avtHistogramFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Allows the filter to change its output's data object information, which
//      is a description of the data object.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Jun 26 09:04:54 PDT 2003
//
//  Modifications:
//    Kathleen Bonnell, Thu Aug 12 08:44:36 PDT 2004
//    Allow this plot to be queryable.
//
// ****************************************************************************

void
avtHistogramFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &outAtts     = GetOutput()->GetInfo().GetAttributes();
    avtDataValidity   &outValidity = GetOutput()->GetInfo().GetValidity();
 
    outAtts.SetTopologicalDimension(1);
    outAtts.SetSpatialDimension(2);
    outValidity.InvalidateZones();
    outValidity.SetNormalsAreInappropriate(true);
    outValidity.InvalidateSpatialMetaData();
    outValidity.SetPointsWereTransformed(true);
}


// ****************************************************************************
//  Method: avtHistogramFilter::ModifyContract
//
//  Purpose:
//      Tells the input that we don't want ghost zones.  They cost time to
//      calculate (in some cases) and skew results towards the domain
//      boundaries.
//
//  Programmer: Hank Childs
//  Creation:   October 18, 2003
//
//  Modifications:
//
//    Hank Childs, Wed Aug 11 09:01:27 PDT 2004
//    Allow the input to have ghost zones, so that we can play better with 
//    other filters.  The histogram filter will now remove ghost zones before
//    executing.
//
//    Hank Childs, Wed May 24 09:44:51 PDT 2006
//    Better support for array variables.
//
//    Hank Childs, Tue Dec 11 20:25:22 PST 2007
//    Add support for weighting by a variable.
//
//    Hank Childs, Tue Feb 19 19:45:43 PST 2008
//    Rename "dynamic" to "streaming", since we really care about whether we
//    are streaming, not about whether we are doing dynamic load balancing.
//    And the two are no longer synonymous.
//
// ****************************************************************************

avtContract_p
avtHistogramFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p newspec = new avtContract(spec);
    if (atts.GetBasedOn() == HistogramAttributes::ManyZonesForSingleVar)
    {
        newspec->NoStreaming();

        if (atts.GetHistogramType() == HistogramAttributes::Variable)
        {
            if (atts.GetWeightVariable() != "default")
                newspec->GetDataRequest()->AddSecondaryVariable(
                                             atts.GetWeightVariable().c_str());
        }
    }
    else
    {
        newspec->GetDataRequest()->TurnZoneNumbersOn();
        avtSILRestriction_p silr = 
                             newspec->GetDataRequest()->GetRestriction();
        int domain = atts.GetDomain();
        int blockOrigin = GetInput()->GetInfo().GetAttributes().GetBlockOrigin();
        domain -= blockOrigin;

        vector<int> domain_list;
        domain_list.push_back(domain);
        silr->TurnOnAll();
        silr->RestrictDomains(domain_list);
    }

    return newspec;
}


// ****************************************************************************
//  Method: avtHistogramFilter::ComputeBinIndex
//
//  Purpose:
//      Compute the bin index for the given input value.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications:
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

int
avtHistogramFilter::ComputeBinIndex( const double &value ) const
{
    // Value out of bounds, just return.
    if ( value < workingRange[0] || value > workingRange[1] )
        return -1;

    // If we have a zero range, return 0
    if ( workingRange[0] == workingRange[1] )
        return 0;
    
    int index = 0;
    if ( atts.GetDataScale() == HistogramAttributes::Linear )
    {
        index = (int)((value - workingRange[0]) / binStep);
    }
    else if ( atts.GetDataScale() == HistogramAttributes::SquareRoot )
    {
        double sign = (value < 0 ? -1.0 : 1.0);
        double x = sign * sqrt( fabs(value) );
        index = (int)((x - sqrtWorkingRange[0]) / sqrtBinStep);        
    }
    else if ( atts.GetDataScale() == HistogramAttributes::Log )
    {
        double sign = (value < 0 ? -1.0 : 1.0);
        double x = sign * log10( fabs( value ) + 1.0 );
        index = (int)((x - logWorkingRange[0]) / logBinStep);
    }

    return index;
}

// ****************************************************************************
//  Method: avtHistogramFilter::ScaleBins
//
//  Purpose:
//      Scale the bins based on the attributes.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications
//    Kevin Bensema Mon Nov 11, 2013. 12:43
//    Added code to compute CDF if desired and to normalize histogram
//    if specified by the attributes.
//
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::ScaleBins()
{
    if ( atts.GetBinScale() == HistogramAttributes::Log )
    {
        for ( int i = 0; i < workingNumBins; i++ )
        {
            double x = bins[i];
            if ( x > 0.0 )
                x = log10(x);
            bins[i] = x;
        }
    }
    else if ( atts.GetBinScale() == HistogramAttributes::SquareRoot )
    {
        for ( int i = 0; i < workingNumBins; i++ )
            bins[i] = sqrt( bins[i] );
    }

  // Normalize the histogram if desired.
  if( atts.GetNormalizeHistogram() == true)
  {
    double sum = 0.0;
    for(int i = 0; i < workingNumBins; ++i)
    {
      sum += bins[i];
    }
    for(int i = 0; i < workingNumBins; ++i)
    {
      bins[i] /= sum;
    }
  }

  // compute the cumulative distribution function if specified in the attributes
  if( atts.GetComputeAsCDF() == true)
  {
    for(int i = 1; i < workingNumBins; ++i)
    {
      bins[i] += bins[i - 1];
    }
  }

  return;
}

// ****************************************************************************
//  Method: avtHistogramFilter::SetWorkingMin
//
//  Purpose:
//      Set the working mins.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications:
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::SetWorkingMin( double dataMin )
{
    workingRange[0] = dataMin;
    double abs_dataMin = fabs( dataMin );
    double sign = (dataMin < 0 ? -1.0 : 1.0);
    logWorkingRange[0] = sign * log10( abs_dataMin + 1.0 );
    sqrtWorkingRange[0] = sign * sqrt( abs_dataMin );
}

// ****************************************************************************
//  Method: avtHistogramFilter::GetWorkingMin
//
//  Purpose:
//      Get the working min based on the type of data scaling.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications:
//    Kathleen Bonnell, Tue Jan  8 17:49:47 PST 2008
//    Added default return.
//
// ****************************************************************************

double
avtHistogramFilter::GetWorkingMin() const
{
    switch ( atts.GetDataScale() )
    {
    case  HistogramAttributes::Linear:
        return workingRange[0];
    case  HistogramAttributes::Log:
        return logWorkingRange[0];
    case HistogramAttributes::SquareRoot:
        return sqrtWorkingRange[0];        
    default:
        return workingRange[0];
    }
}

// ****************************************************************************
//  Method: avtHistogramFilter::SetWorkingMax
//
//  Purpose:
//      Set the working maxs.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications:
//    Mark C. Miller, Thu Jan  8 14:27:03 PST 2015
//    Adjusted to use doubles everywhere.
// ****************************************************************************

void
avtHistogramFilter::SetWorkingMax( double dataMax )
{
    workingRange[1] = dataMax;
    double abs_dataMax = fabs( dataMax );

    double sign = (dataMax < 0 ? -1.0 : 1.0);
    logWorkingRange[1] = sign * log10( abs_dataMax + 1.0 );
    sqrtWorkingRange[1] = sign * sqrt( abs_dataMax );
}

// ****************************************************************************
//  Method: avtHistogramFilter::GetWorkingMax
//
//  Purpose:
//      Get the working max based on the type of data scaling.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
//  Modifications:
//    Kathleen Bonnell, Tue Jan  8 17:49:47 PST 2008
//    Added default return.
//
// ****************************************************************************

double
avtHistogramFilter::GetWorkingMax() const
{
    switch ( atts.GetDataScale() )
    {
    case  HistogramAttributes::Linear:
        return workingRange[1];
    case  HistogramAttributes::Log:
        return logWorkingRange[1];
    case HistogramAttributes::SquareRoot:        
        return sqrtWorkingRange[1];        
    default :
        return workingRange[1];
    }
}

// ****************************************************************************
//  Method: avtHistogramFilter::SetWorkingNumBins
//
//  Purpose:
//      Set the number of working bins and the bin step sizes.
//
//  Programmer: Dave Pugmire
//  Creation:   November 01, 2007
//
// ****************************************************************************

void
avtHistogramFilter::SetWorkingNumBins( int n )
{
    workingNumBins = n;
    binStep = (workingRange[1] - workingRange[0]) / workingNumBins;
    logBinStep = (logWorkingRange[1] - logWorkingRange[0] ) / workingNumBins;
    sqrtBinStep = (sqrtWorkingRange[1] - sqrtWorkingRange[0] ) / workingNumBins;
}
